#Load data

Field data, with predictions with INLA; interpolations with INLA
data from IBM model

#packages
library(tidyr)
library(dplyr)
library(INLA)
library(plotly)
library(forecast)
library(gridExtra)

#field data and predictions from INLA
lidar_data <- read.csv("output/Data_and_predictions_INLA_fielddata.csv", sep=",")
interlidar_data <- read.csv("output/interpolations_INLA_fielddata.csv", sep=",")

#data IBM output
ibm_data <- read.csv("data/data_model.csv", sep=",") %>% slice(-836) %>% #delete wrong (double) line
  filter(error < 0.4) %>% #filter outlier from error
  mutate(rep = as.factor(rep))

#Transform variables

To get comparable values/units and scale them the same

#summary(lidar_data)
#summary(ibm_data)

#divide JC in two, than comparable to field data
ibm_data <- ibm_data %>% mutate(JC2 = JC/2) %>% 
  #standardize IBM data: P is from 0 to 1, so ok; Hstd scaled JC2
  mutate(Pstd = P, Hstd = scale(JC2, center=mean(JC2), scale=sd(JC2)))

#Models for IBM data

#first, get data for interpolations
interibm_data <- tibble(Pstd=seq(0,1, by=0.01), H=seq(20,70,by=0.5)) %>%
  tidyr::expand(Pstd, H) %>%
  mutate(Hstd=scale(H, center=mean(ibm_data$JC2), scale=sd(ibm_data$JC2))) %>% #scale by IBM data, to get similar values
  mutate(volume = NA, error = NA)#add empty values for volume and error
  
#select needed columns from ibm_data
ibm_all <- ibm_data %>% select(Pstd, Hstd, volume, error) %>%
  #bind with interpol data to feed it to INLA
  bind_rows(select(interibm_data, -H))
  
#INLA
#for volume
f0_diff <- volume ~  Pstd + Hstd + Pstd:Hstd +
  I(Pstd^2) + I(Hstd^2) +Hstd:I(Pstd^2) + Pstd:I(Hstd^2)
M0_diff <- inla(f0_diff, 
                #control.compute = list(dic = TRUE, waic = TRUE),
                control.predictor = list(compute = TRUE),
                family = "gaussian", data = ibm_all)
summary(M0_diff)
## 
## Call:
## c("inla(formula = f0_diff, family = \"gaussian\", data = ibm_all, ",  "    control.predictor = list(compute = TRUE))")
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          0.6524         10.7607          1.1690         12.5822 
## 
## Fixed effects:
##                     mean      sd 0.025quant  0.5quant 0.975quant      mode
## (Intercept)      55.8184  6.6912    42.7096   55.8082    68.9719   55.7884
## Pstd             84.3977 20.4905    44.0965   84.4218   124.5331   84.4720
## Hstd            -31.6104  5.7107   -42.8082  -31.6158   -20.3933  -31.6261
## I(Pstd^2)      -132.8997 19.4652  -171.0398 -132.9258   -94.6450 -132.9768
## I(Hstd^2)        -4.3092  4.6528   -13.4432   -4.3102     4.8221   -4.3118
## Pstd:Hstd        68.0558 19.5353    29.6698   68.0658   106.3521   68.0877
## Hstd:I(Pstd^2)  -12.2823 19.5982   -50.7301  -12.2932    26.1917  -12.3136
## Pstd:I(Hstd^2)    9.5700  7.8981    -5.9401    9.5700    25.0652    9.5708
##                kld
## (Intercept)      0
## Pstd             0
## Hstd             0
## I(Pstd^2)        0
## I(Hstd^2)        0
## Pstd:Hstd        0
## Hstd:I(Pstd^2)   0
## Pstd:I(Hstd^2)   0
## 
## The model has no random effects
## 
## Model hyperparameters:
##                                          mean   sd 0.025quant 0.5quant
## Precision for the Gaussian observations 3e-04 0.00      2e-04    3e-04
##                                         0.975quant  mode
## Precision for the Gaussian observations      3e-04 3e-04
## 
## Expected number of effective parameters(std dev): 6.337(0.0272)
## Number of equivalent replicates : 156.06 
## 
## Marginal log-Likelihood:  -5541.16 
## Posterior marginals for linear predictor and fitted values computed
#for error
f0_err <- error ~ Pstd + Hstd + Pstd:Hstd +
  I(Pstd^2) + I(Hstd^2) +Hstd:I(Pstd^2) + Pstd:I(Hstd^2)
M0_err <- inla(f0_err,
               control.predictor = list(compute = TRUE),
               family = "gamma", data = ibm_all)
summary(M0_err)
## 
## Call:
## "inla(formula = f0_err, family = \"gamma\", data = ibm_all, control.predictor = list(compute = TRUE))"
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          0.5360         61.3286          1.4010         63.2656 
## 
## Fixed effects:
##                   mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept)    -4.5497 0.0569    -4.6608  -4.5499    -4.4375 -4.5504   0
## Pstd            0.9001 0.2278     0.4518   0.9005     1.3460  0.9012   0
## Hstd           -0.1399 0.0511    -0.2401  -0.1399    -0.0395 -0.1400   0
## I(Pstd^2)      -1.2839 0.2134    -1.7022  -1.2842    -0.8644 -1.2847   0
## I(Hstd^2)       0.0284 0.0304    -0.0311   0.0283     0.0881  0.0282   0
## Pstd:Hstd       0.9552 0.2226     0.5177   0.9553     1.3915  0.9556   0
## Hstd:I(Pstd^2) -0.6835 0.2213    -1.1178  -0.6837    -0.2489 -0.6839   0
## Pstd:I(Hstd^2)  0.0038 0.0520    -0.0983   0.0038     0.1057  0.0038   0
## 
## The model has no random effects
## 
## Model hyperparameters:
##                                                 mean     sd 0.025quant
## Precision parameter for the Gamma observations 7.021 0.3052      6.435
##                                                0.5quant 0.975quant  mode
## Precision parameter for the Gamma observations    7.016      7.634 7.007
## 
## Expected number of effective parameters(std dev): 8.042(0.0019)
## Number of equivalent replicates : 122.98 
## 
## Marginal log-Likelihood:  3981.09 
## Posterior marginals for linear predictor and fitted values computed
save(M0_diff, file="output/M0_diff.rda")
save(M0_err, file="output/M0_err.rda")

#extract interpolation data from IBM models
M0_diff_inter <- M0_diff$summary.fitted.values %>%
  slice(990:n()) %>% #select predictions from dataset that are not input data
  select(mean, sd, `0.025quant`, `0.975quant`) %>% #select neaded values
  #rename to get similar values as in interlidar_data
  rename(Int_diff_mean = mean, Int_diff_sd = sd, Int_diff_0.025quant=`0.025quant`, Int_diff_0.975quant = `0.975quant`)

M0_err_inter <- M0_err$summary.fitted.values %>%
  slice(990:n()) %>% #select predictions from dataset that are not input data
  select(mean, sd, `0.025quant`, `0.975quant`) %>% #select neaded values
  #rename to get similar values as in interlidar_data
  rename(Int_err_mean = mean, Int_err_sd = sd, Int_err_0.025quant=`0.025quant`, Int_err_0.975quant = `0.975quant`)



#bind these interpolations/predictions to interibm_data
interibm_data <- interibm_data %>% bind_cols(M0_diff_inter) %>% bind_cols(M0_err_inter)

#Effect sizes

Visualisation of the effect sizes/signs of the covariates

#load INLA models from lidar data
load("output/M1_diff.rda")
load("output/M1_err.rda")

#lidar data get coefficients from models
post_beta1_diff <- M1_diff$summary.fixed[, c("mean", "0.025quant", "0.975quant")]
post_beta1_err <- M1_err$summary.fixed[, c("mean", "0.025quant", "0.975quant")]

post_beta1_diff <- post_beta1_diff %>% tibble::rownames_to_column(var="WhichVar") %>%
  mutate(WhichVar = factor(WhichVar, levels=c("meandVstd", "Diststd",  "Pstd:I(Hstd^2)", "Hstd:I(Pstd^2)",
                                              "I(Hstd^2)","I(Pstd^2)", "Pstd:Hstd", "Hstd", "Pstd","Intercept" ))) %>%
  mutate(WhichVar = recode(WhichVar, "meandVstd"="meandV", "Diststd"="Dist",
                           "Pstd:I(Hstd^2)"="P:JC²", "Hstd:I(Pstd^2)"="JC:P²",
                           "I(Hstd^2)"="JC²","I(Pstd^2)"="P²", "Pstd:Hstd"="P:JC", "Hstd"="JC",
                           "Pstd"="P","(Intercept)"="(Intercept)" ))

post_beta1_err <- post_beta1_err %>% tibble::rownames_to_column(var="WhichVar") %>%
  mutate(WhichVar = factor(WhichVar, levels=c("meandVstd", "Diststd",  "Pstd:I(Hstd^2)", "Hstd:I(Pstd^2)",
                                              "I(Hstd^2)","I(Pstd^2)", "Pstd:Hstd", "Hstd", "Pstd","Intercept" )))  %>%
  mutate(WhichVar = recode(WhichVar, "meandVstd"="meandV", "Diststd"="Dist",
                           "Pstd:I(Hstd^2)"="P:JC²", "Hstd:I(Pstd^2)"="JC:P²",
                           "I(Hstd^2)"="JC²","I(Pstd^2)"="P²", "Pstd:Hstd"="P:JC", "Hstd"="JC",
                           "Pstd"="P","(Intercept)"="(Intercept)" ))

#ibm data get coefficients from models
post_beta0_diff <- M0_diff$summary.fixed[, c("mean", "0.025quant", "0.975quant")]
post_beta0_err <- M0_err$summary.fixed[, c("mean", "0.025quant", "0.975quant")]

post_beta0_diff <- post_beta0_diff %>% tibble::rownames_to_column(var="WhichVar") %>%
  mutate(WhichVar = factor(WhichVar, levels=c("Pstd:I(Hstd^2)", "Hstd:I(Pstd^2)",
                                              "I(Hstd^2)","I(Pstd^2)", "Pstd:Hstd", "Hstd", "Pstd","(Intercept)" ))) %>%
  mutate(WhichVar = recode(WhichVar, "Pstd:I(Hstd^2)"="P:JC²", "Hstd:I(Pstd^2)"="JC:P²",
         "I(Hstd^2)"="JC²","I(Pstd^2)"="P²", "Pstd:Hstd"="P:JC", "Hstd"="JC",
         "Pstd"="P","(Intercept)"="(Intercept)" ))

post_beta0_err <- post_beta0_err %>% tibble::rownames_to_column(var="WhichVar") %>%
  mutate(WhichVar = factor(WhichVar, levels=c("Pstd:I(Hstd^2)", "Hstd:I(Pstd^2)",
                                              "I(Hstd^2)","I(Pstd^2)", "Pstd:Hstd", "Hstd", "Pstd","(Intercept)" ))) %>%
  mutate(WhichVar = recode(WhichVar, "Pstd:I(Hstd^2)"="P:JC²", "Hstd:I(Pstd^2)"="JC:P²",
         "I(Hstd^2)"="JC²","I(Pstd^2)"="P²", "Pstd:Hstd"="P:JC", "Hstd"="JC",
         "Pstd"="P","(Intercept)"="(Intercept)" ))

#lidar data
p_diff_lidar <- ggplot(post_beta1_diff) +
  geom_pointrange() +
  geom_hline(yintercept = 0, linetype = 2) +
  aes(x=WhichVar, y =mean, ymin = `0.025quant`, ymax = `0.975quant`) +
  xlab("covariate")+ ylab("Effect size")+
  ggtitle("Posteriors for \u0394height (LiDAR data)")+
  theme_classic() +
  coord_flip()
#p_diff_lidar

p_err_lidar <- ggplot(post_beta1_err) +
  geom_pointrange() +
  geom_hline(yintercept = 0, linetype = 2) +
  aes(x=WhichVar, y =mean, ymin = `0.025quant`, ymax = `0.975quant`) +
  xlab("covariate")+ ylab("Effect size")+
  ggtitle("Posteriors for CV(\u0394height) (LiDAR data)")+
  theme_classic() +
  coord_flip()
#p_err_lidar

#IBM model
p_diff_ibm <- ggplot(post_beta0_diff) +
  geom_pointrange() +
  geom_hline(yintercept = 0, linetype = 2) +
  aes(x=WhichVar, y =mean, ymin = `0.025quant`, ymax = `0.975quant`) +
  xlab("covariate")+ ylab("Effect size")+
  ggtitle("Posteriors for volume (IBM data)")+
  theme_classic() +
  coord_flip()
#p_diff_ibm

p_err_ibm <- ggplot(post_beta0_err) +
  geom_pointrange() +
  geom_hline(yintercept = 0, linetype = 2) +
  aes(x=WhichVar, y =mean, ymin = `0.025quant`, ymax = `0.975quant`) +
  xlab("covariate")+ ylab("Effect size")+
  ggtitle("Posteriors for CV(volume) (IBM data)")+
  theme_classic() +
  coord_flip()
#p_err_ibm

grid.arrange(p_diff_ibm, p_err_ibm, p_diff_lidar, p_err_lidar)

#Visualisation of the interpolations

3D graphs for both the model and lidar data

Difference (height/volume): lijkt me bij combinate hoge P en H en kleine P en H enkel verschillend van elkaar. Bij lidar data gaat bij hoge P en H gaat difference ferm naar beneden, bij IBM data niet. Bij lage P en H gaat bij ibm data de difference naar omhoog, bij lidar data weer naar beneden. Bij de niet extreme waarden lijkt het wel goed overeen te komen (zou eens apart kunnen geplot worden).

Error: zadelprofiel bij ibm data, bij lidar data niet aanwezig (grote en kleine waarden van H)…

#3Dplot
#Lidar data
Interlidar_diff_3D <- plot_ly(data = interlidar_data, x = ~H, y = ~Pstd, z=~Int_diff_mean, color=~Int_diff_mean)  %>% add_markers() %>%
  layout(title = "Difference in dune height (LiDAR data)",
  scene = list(
    xaxis = list(title = "JC"),
    yaxis = list(title = "P", tickvals=c(0,0.2,0.4,0.6,0.8,1)),
    zaxis = list(title = "\u0394height (m)", titleangle=90))) %>%
  colorbar(title = "Prediction \u0394height")
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
Interlidar_err_3D <- plot_ly(data = interlidar_data, x = ~H, y = ~Pstd, z=~exp(Int_err_mean), color=~exp(Int_err_mean))  %>% add_markers() %>%
  layout(title = "Variability in dune height (LiDAR data)",
  scene = list(
    xaxis = list(title = "JC"),
    yaxis = list(title = "P", tickvals=c(0,0.2,0.4,0.6,0.8,1)),
    zaxis = list(title = "CV (m)"))) %>%
  colorbar(title = "Prediction CV")



#IBM model
Interibm_diff_3D <- plot_ly(data = interibm_data, x = ~H, y = ~Pstd, z=~Int_diff_mean, color=~Int_diff_mean)  %>% add_markers() %>%
  layout(title = "Difference in dune volume (IBM data)",
  scene = list(
    xaxis = list(title = "JC"),
    yaxis = list(title = "P", tickvals=c(0,0.2,0.4,0.6,0.8,1)),
    zaxis = list(title = "volume (m³)"))) %>%
  colorbar(title = "Prediction volume")


Interibm_err_3D <- plot_ly(data = interibm_data, x = ~H, y = ~Pstd, z=~exp(Int_err_mean), color=~exp(Int_err_mean))  %>% add_markers() %>%
  layout(title = "Variability in dune volume (IBM data)",
  scene = list(
    xaxis = list(title = "JC"),
    yaxis = list(title = "P", tickvals=c(0,0.2,0.4,0.6,0.8,1)),
    zaxis = list(title = "CV (m³)"))) %>%
  colorbar(title = "Prediction CV")



Interlidar_diff_3D
Interibm_diff_3D
Interlidar_err_3D
Interibm_err_3D

#Plotting uncertainties

#3Dplot
#Lidar data
Interlidar_diff_3D_uncertainty <- plot_ly(data = interlidar_data, x = ~H, y = ~Pstd, z=~Int_diff_mean, color=~(Int_diff_0.975quant - Int_diff_0.025quant))  %>% add_markers() %>%
  layout(title = "Uncertainty Difference (lidar data)",
  scene = list(
    xaxis = list(title = "JC"),
    yaxis = list(title = "P", tickvals=c(0,0.2,0.4,0.6,0.8,1)),
    zaxis = list(title = "Height (m)")))

Interlidar_err_3D_uncertainty <- plot_ly(data = interlidar_data, x = ~H, y = ~Pstd, z=~exp(Int_err_mean), color=~exp(Int_err_0.975quant-Int_err_0.025quant))  %>% add_markers() %>%
  layout(title = "Uncertainty Variability (lidar data)",
  scene = list(
    xaxis = list(title = "JC"),
    yaxis = list(title = "P"),
    zaxis = list(title = "Error (m)")))

Interibm_diff_3D_uncertainty <- plot_ly(data = interibm_data, x = ~H, y = ~Pstd, z=~Int_diff_mean, color=~(Int_diff_0.975quant - Int_diff_0.025quant))  %>% add_markers() %>%
  layout(title = "Uncertainty volume (ibm data)",
  scene = list(
    xaxis = list(title = "JC"),
    yaxis = list(title = "P"),
    zaxis = list(title = "Height (m)")))

Interibm_err_3D_uncertainty <- plot_ly(data = interibm_data, x = ~H, y = ~Pstd, z=~exp(Int_err_mean), color=~exp(Int_err_0.975quant-Int_err_0.025quant))  %>% add_markers() %>%
  layout(title = "Uncertainty variability (ibm data)",
  scene = list(
    xaxis = list(title = "JC"),
    yaxis = list(title = "P"),
    zaxis = list(title = "Error (m)")))


Interlidar_diff_3D_uncertainty
Interibm_diff_3D_uncertainty
Interlidar_err_3D_uncertainty
Interibm_err_3D_uncertainty